First, let’s look at raster metadata using GDALinfo

GDALinfo("./data/HARV_dsmCrop.tif")

rows        1367 
columns     1697 
bands       1 
lower left origin.x        731453 
lower left origin.y        4712471 
res.x       1 
res.y       1 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
file        ./data/HARV_dsmCrop.tif 
apparent band summary:
   GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1 Float64           TRUE       -9999          1       1697
apparent band statistics:
    Bmin   Bmax    Bmean      Bsd
1 305.07 416.07 359.8531 17.83169
Metadata:
AREA_OR_POINT=Area 

Load raster

Note: raster() will create a single-band raster. stack() and brick() can be used to load/create a multi-band raster (more on that later)

harv_dsm <- raster("./data/HARV_dsmCrop.tif")
harv_dsm

class       : RasterLayer 
dimensions  : 1367, 1697, 2319799  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 731453, 733150, 4712471, 4713838  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /Users/francesdavenport/Documents/teaching/r_geospatial/r-geospatial-102219/data/HARV_dsmCrop.tif 
names       : HARV_dsmCrop 
values      : 305.07, 416.07  (min, max)

We can plot a raster using base R plotting…

plot(harv_dsm)

…or using ggplot (after creating a data.frame with the data from our raster)

harv_dsm_df <- as.data.frame(harv_dsm, xy = TRUE)
str(harv_dsm_df)

'data.frame':   2319799 obs. of  3 variables:
 $ x           : num  731454 731454 731456 731456 731458 ...
 $ y           : num  4713838 4713838 4713838 4713838 4713838 ...
 $ HARV_dsmCrop: num  409 408 407 407 409 ...

ggplot() + 
  geom_raster(data = harv_dsm_df, aes(x, y, fill = HARV_dsmCrop)) + 
  scale_fill_viridis_c()

Let’s calculate canopy height from a DSM (digital surface model) and a DTM (digital terrain model)

Let’s read in our digital terrain model

harv_dtm <- raster("./data/HARV_dtmCrop.tif")
plot(harv_dtm)

Raster math: we can perform calculations between two rasters

Let’s calculate canopy height by subtracting the DTM from the DSM

harv_canopy <- harv_dsm - harv_dtm
plot(harv_canopy)

Raster math can be slow for large rasters or complex calculations…an alternative is to use overlay

harv_canopy <- overlay(harv_dsm, harv_dtm, 
                       fun = function(r1, r2) { return( r1 - r2) })
plot(harv_canopy)

We can also do calculations between rasters and integers

Let’s calculate the mean height in our DTM and subtract that from our DTM

harv_dtm2 <- harv_dtm - cellStats(harv_dtm, stat='mean')
plot(harv_dtm2)

Raster calc with logical expressions

harv_dtm3 <- calc(harv_dtm, function(x) ifelse(x > 340, x, NA))
plot(harv_dtm3)

Saving our canopy height model for future use:

Let’s write a GeoTiff file. R will recognize the desired output file type based on the file extension (see documentation for options)

writeRaster(harv_canopy, "./data/HARV_canopy_height.tif")

Exercise: come up with your own raster calculation and see if it works! Discuss with your neighbor

Let’s read in a shapefile with the boundary of our field site:

Using the sp package:

aoi_sp <- shapefile("./data/HarClip_UTMZ18.shp")
aoi_sp

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 732128, 732251.1, 4713209, 4713359  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 1
names       : id 
value       :  1 

plot(aoi_sp)

Now let’s read in the same shapefile using the sf (“simple features”) package:

library(sf)

Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

aoi <- st_read("./data/HarClip_UTMZ18.shp")

Reading layer `HarClip_UTMZ18' from data source `/Users/francesdavenport/Documents/teaching/r_geospatial/r-geospatial-102219/data/HarClip_UTMZ18.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 732128 ymin: 4713209 xmax: 732251.1 ymax: 4713359
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

plot(aoi)

Let’s try another shapefile with multiple objects:

ne_states <- st_read("./data/Boundary-US-State-NEast.shp")

Reading layer `Boundary-US-State-NEast' from data source `/Users/francesdavenport/Documents/teaching/r_geospatial/r-geospatial-102219/data/Boundary-US-State-NEast.shp' using driver `ESRI Shapefile'
Simple feature collection with 12 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XYZ
bbox:           xmin: -80.51989 ymin: 37.91685 xmax: -66.9499 ymax: 47.45716
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

head(ne_states)

Simple feature collection with 6 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XYZ
bbox:           xmin: -79.76212 ymin: 37.91685 xmax: -66.9499 ymax: 47.45716
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  STATEFP  STATENS    AFFGEOID GEOID STUSPS                 NAME LSAD        ALAND      AWATER
1      11 01702382 0400000US11    11     DC District of Columbia   00    158350578    18633500
2      24 01714934 0400000US24    24     MD             Maryland   00  25147575220  6983455225
3      36 01779796 0400000US36    36     NY             New York   00 122054577774 19242052501
4      09 01779780 0400000US09    09     CT          Connecticut   00  12542396439  1814978794
5      23 01779787 0400000US23    23     ME                Maine   00  79885244456 11749091526
6      25 00606926 0400000US25    25     MA        Massachusetts   00  20203936287  7131805598
                        geometry
1 MULTIPOLYGON Z (((-77.11976...
2 MULTIPOLYGON Z (((-76.04621...
3 MULTIPOLYGON Z (((-72.01893...
4 MULTIPOLYGON Z (((-73.69594...
5 MULTIPOLYGON Z (((-68.92401...
6 MULTIPOLYGON Z (((-70.27553...

plot(ne_states)

names(ne_states)

 [1] "STATEFP"  "STATENS"  "AFFGEOID" "GEOID"    "STUSPS"   "NAME"     "LSAD"     "ALAND"    "AWATER"   "geometry"

We can plot SF objects with ggplot!

ggplot() + 
  geom_sf(data = ne_states, aes(fill = NAME))

Some common sf-specific commands:

st_bbox(ne_states)

     xmin      ymin      xmax      ymax 
-80.51989  37.91685 -66.94989  47.45716 

st_crs(ne_states)

Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Simple features objects also work with many dplyr/tidyr type commands:

nrow(ne_states)

[1] 12

mass <- ne_states %>%
  filter(NAME == "Massachusetts")
ggplot() + 
  geom_sf(data = mass, aes(fill = NAME))

In addition to polygon shapefiles, we can read point and line data:

tower <- st_read("./data/HARVtower_UTM18N.shp")

Reading layer `HARVtower_UTM18N' from data source `/Users/francesdavenport/Documents/teaching/r_geospatial/r-geospatial-102219/data/HARVtower_UTM18N.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 14 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 732183.2 ymin: 4713265 xmax: 732183.2 ymax: 4713265
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

roads <- st_read("./data/HARV_roads.shp")

Reading layer `HARV_roads' from data source `/Users/francesdavenport/Documents/teaching/r_geospatial/r-geospatial-102219/data/HARV_roads.shp' using driver `ESRI Shapefile'
Simple feature collection with 13 features and 15 fields
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: 730741.2 ymin: 4711942 xmax: 733295.5 ymax: 4714260
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

site_map <- ggplot(roads) + 
  geom_sf() + 
  geom_sf(data = tower, col = "red", size = 10, shape = "*")

Note: ggplot is displaying coordinates as lat/lon, even though our vector data has the UTM projection!

st_crs(tower)

Coordinate Reference System:
  EPSG: 32618 
  proj4string: "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs"

Create 200m buffer around flux tower

tower_buffer <- st_buffer(tower, 200)
site_map + 
  geom_sf(data = tower_buffer, fill = "green", alpha = 0.2) 

“Clip” canopy height model to 200m around flux tower

Note: We could use crop to clip the raster to the bounding box of an object, but not it’s actual boundary

canopy_masked <- mask(harv_canopy, tower_buffer)
canopy_masked

class       : RasterLayer 
dimensions  : 1367, 1697, 2319799  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 731453, 733150, 4712471, 4713838  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 38.16998  (min, max)

plot(canopy_masked)

Instead, maybe we just want the canopy height values for all cells within our buffer

buffer_canopy_values <- extract(harv_canopy, tower_buffer)

Extract returns a list of vectors (one for each polygon object), but we only have one object so we can unlist the result

buffer_canopy_values <- unlist(buffer_canopy_values)
hist(buffer_canopy_values)

Let’s find our NDVI (normalized difference vegetation index) files and read them in

ndvi_files <- list.files("./data/NDVI/", 
                            full.names = TRUE,
                            pattern = ".tif$")
harv_ndvi <- stack(ndvi_files)
harv_ndvi

class       : RasterStack 
dimensions  : 5, 4, 20, 13  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent      : 239415, 239535, 4714215, 4714365  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=19 +ellps=WGS84 +units=m +no_defs 
names       : X005_HARV_ndvi_crop, X037_HARV_ndvi_crop, X085_HARV_ndvi_crop, X133_HARV_ndvi_crop, X181_HARV_ndvi_crop, X197_HARV_ndvi_crop, X213_HARV_ndvi_crop, X229_HARV_ndvi_crop, X245_HARV_ndvi_crop, X261_HARV_ndvi_crop, X277_HARV_ndvi_crop, X293_HARV_ndvi_crop, X309_HARV_ndvi_crop 
min values  :                2732,                1534,                1917,                5669,                8685,                8768,                8633,                8686,                8297,                7491,                 277,                 484,                4396 
max values  :                5545,                3736,                3600,                6394,                8903,                9140,                8945,                8967,                8651,                8607,                 366,                 661,                6359 

crs(harv_ndvi)

CRS arguments:
 +proj=utm +zone=19 +ellps=WGS84 +units=m +no_defs 

Wait, this is a different zone! Let’s reproject the raster to UTM zone 18

crs(harv_dtm)

CRS arguments:
 +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 

harv_ndvi <- projectRaster(harv_ndvi, crs = crs(harv_dtm))

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

Now, let’s calculate mean NDVI on each date

ndvi_mean <- cellStats(harv_ndvi, mean)
plot(ndvi_mean)

Instead, let’s calculate mean NDVI at our tower location

We could do this by finding the cell overlapping with the tower and extracting the value at that point

tower_xy <- st_coordinates(tower)
tower_row <- rowFromY(harv_ndvi, tower_xy[2])
tower_col <- colFromX(harv_ndvi, tower_xy[1])
tower_ndvi <- getValues(harv_ndvi, row = tower_row)[tower_col,]
plot(tower_ndvi)

We can do the same thing using extract (with fewer lines of code!)

tower_ndvi <- raster::extract(harv_ndvi, tower)
buffer_ndvi <- raster::extract(harv_ndvi, tower_buffer, fun = mean)
buffer_ndvi

     X005_HARV_ndvi_crop X037_HARV_ndvi_crop X085_HARV_ndvi_crop X133_HARV_ndvi_crop X181_HARV_ndvi_crop
[1,]            3709.251            2507.582            2571.568            6004.585             8780.76
     X197_HARV_ndvi_crop X213_HARV_ndvi_crop X229_HARV_ndvi_crop X245_HARV_ndvi_crop X261_HARV_ndvi_crop
[1,]            8941.279            8783.959            8818.749            8503.574            7955.772
     X277_HARV_ndvi_crop X293_HARV_ndvi_crop X309_HARV_ndvi_crop
[1,]            332.0345            561.2631            5457.627

Optional: Let’s make this data look a little bit nicer and plot it:

ndvi_ts <- as.data.frame(unlist(buffer_ndvi)[1,])
names(ndvi_ts) <- "meanNDVI"
julian_day <- gsub("X|_HARV_ndvi_crop", "", row.names(ndvi_ts))
origin <- as.Date("2011-01-01")
ndvi_ts$date <- origin + as.integer(julian_day) - 1
ndvi_ts

                     meanNDVI       date
X005_HARV_ndvi_crop 3709.2514 2011-01-05
X037_HARV_ndvi_crop 2507.5816 2011-02-06
X085_HARV_ndvi_crop 2571.5682 2011-03-26
X133_HARV_ndvi_crop 6004.5853 2011-05-13
X181_HARV_ndvi_crop 8780.7595 2011-06-30
X197_HARV_ndvi_crop 8941.2789 2011-07-16
X213_HARV_ndvi_crop 8783.9587 2011-08-01
X229_HARV_ndvi_crop 8818.7488 2011-08-17
X245_HARV_ndvi_crop 8503.5736 2011-09-02
X261_HARV_ndvi_crop 7955.7723 2011-09-18
X277_HARV_ndvi_crop  332.0345 2011-10-04
X293_HARV_ndvi_crop  561.2631 2011-10-20
X309_HARV_ndvi_crop 5457.6267 2011-11-05

ggplot(ndvi_ts, aes(date, meanNDVI)) + 
  geom_point() + 
  geom_line()

Exercise: Come up with your own analysis using the raster and shapefile data that we have

Remember, our vector data includes roads and our area of interest, in addition to making our own areas using things like st_buffer! Try to plot your results and discuss it with your neighbor